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, Q ' Abstract 

. A new DRP scheme is built, which enables us to minimize the error due to 

OO . 

, the finite difference approximation, by means of an equivalent matrix equation. 
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2 ; 1 Introduction: Scheme classes 



We hereafter propose a method that enables us to build a DRP scheme while mini- 
^ , mizing the error due to the finite difference approximation, by means of an equivalent 

\^ ' matrix equation. 
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Consider the transport equation: 



■ dU du r n r 

^ \ with the initial condition u{x,t = 0) = Uq{x). 
> 

^ ' Proposition 1.1 A finite difference scheme for this equation can be written under 

■ the form: 

a u,"+7 m^-i^+C u^_^^'^+9 u^+^+d m+i"-' = 

(2) 

where: 

ur =u{lh,mT) (3) 

/ G — i, « + fn E {n — 1, n, n + 1}, j = 0, n^, n = 0, rit, h, r denoting 
respectively the mesh size and time step (L = Ux h, T = ritr). 
The Courant-Friedrichs-Lewy number (cfl) is defined as a = cr/h . 
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Table 1: Numerical scheme coefficient. 
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A numerical scheme is specified by selecting appropriate values of the coefficients 
a, P, 7, 6, e, (, r], 6 and in equation ([2]), which, for sake of usefulness, will be 
written as: 

a = a^^ + at , (3 = (3^^ + (3t , 7 = 7x + 7t , 5 = 5x + 5t , e = e^^ + St , (4) 

where the "a;" denotes a dependance towards the mesh size h, while the "f" denotes 
a dependance towards the time step r. 

Values corresponding to numerical schemes retained for the present works are given 
in Table [H 

The number of time steps will be denoted n^, the number of space steps, Ux- In 
general, ^ rit. 

In the following: the only dependance of the coefficients towards the time step r 
existing only in the Crank-Nicolson scheme, we will restrain our study to the specific 
case: 

«t = 7t = C = ^ = ^ = ^ = (5) 

The paper is organized as follows. The building of the DRP scheme is exposed in 
section [21 The equivalent matrix equation, which enables us to minimize the error 
due to the finite difference approximation, is presented in section [31 A numerical 
example is given in section [H 

2 The DRP scheme 

The first derivative |^ is approximated at the node of the spatial mesh by: 

du 

( ^ )« ^ /5x wj+i" + 5x + Wi+i-i" (6) 

Following the method exposed by C. Tam and J. Webb in [1], the coefficients (3xi 
8x1 and Ex are determined requiring the Fourier Transform of the finite difference 
scheme ([6]) to be a close approximation of the partial derivative ( )«■ 
([6]) is a special case of: 



DRP SCHEME OPTIMIZATION 



3 



du 

{ )ic^p^u{x + ih) + 6x u{x + {i + l)h)+ Ex u{x + - 1) h) (7) 
ox 

where x is a continuous variable, and can be recovered setting x = I h. 
Applying the Fourier transform, referred to by ^, to both sides of ([7]), yields: 

{/3,.e° + 5,.e^^-^ + £,e-^'-^} u (8) 
j denoting the complex square root of —1. 

Comparing the two sides of ([H]) enables us to identify the wavenumber A of the finite 
difference scheme and the quantity j [P^^^ + ^x^'^'^^ + £x€-~^'^^}, i. e.: The 
wavenumber of the finite difference scheme dH]) is thus: 

A = -J {13x6^ + 6x6^''^ + 6x6'^^^} (9) 

To ensure that the Fourier transform of the finite difference scheme is a good ap- 
proximation of the partial derivative ( |^ )/ over the range of waves with wavelength 
longer than 4/i, the a priori unknowns coefficients (3x, Sx, and Ex must be choosen 
so as to minimize the integrated error: 



/A \Xh- Xh\'^d{Xh) 

jK \K + jh {Pxe^ + Sxe^^' + exe-^^} ^^(k) 



(10) 



The conditions that £^ is a minimum are: 

^^K^K^O (11) 

d/3x d6x dsx 

and provide the following system of linear algebraic equations: 

2Txhl3x + 4:{h5x + he^-l) = 

Ah(3x + 7r{2 6x-l) = (12) 
4:hpx + '2Trh€x = 

which enables us to determine the required values of Px, ^x, and Ex'. 
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= 5r 
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^x 



/i(7r2-8) 
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(13) 



2 h{n^-8) 
2 

/i(7r2-8) 

3 The Sylvester equation 

3.1 Matricial form of the finite diff'erences problem 

Theorem 3.1 The problem ([2]) can be written under the following matricial form: 

M1U + UM2 + C{U) = Mo (14) 
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where Mi and M2 are square matrices respectively Ux — 1 hy Ux — I, rit hy rit, given 
by: 



Ml 



f P S ... \ 

e (3 ■-. ■-. : 

■-. ■-. ■-. 

: ■-. ■-. /3 s 

V ... e P / 



the matrix Mq being given by: 



Mo 



/ 7 ... \ 

a ■-. ■-. : 

■-. ■-. ■-. 

: ■-. ■-. ■-. 7 

V ... a J 



(15) 
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(16) 



and where £ is a hnear matricial operator which can be written as: 

£ = £1 + £2 + A + A 
where £1, £2, and £4 are given by: 



cm = c 
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(19) 



Proposition 3.2 The second member matrix Mq bears the initial conditions, given 
for the specific value n — 0, which correspond to the initialization process when 
computing loops, and the boundary conditions, given for the specific values i — 0, 
i = Ux- 
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Denote by Uexact the exact solution of ([T]). 
The corresponding matrix Uexact will be: 

Uexact [Uexacti ] l<i<nx — l,l<n<nt (^0) 

where: 



Uexacti Uexacti^-^i^ ^n) (^-^) 



with Xi = i h, tn = n T. 



Definition 3.3 We will call error matrix the matrix defined by: 

E = U~ Uexact (22) 

Consider the matrix F defined by: 

F = Ml Uexact + Uexact + C{Uexact) ' Mq (23) 

Proposition 3.4 The error matrix E satisfies: 

MiE + EM2 + C{E) = F (24) 

3.2 The matrix equation 

Theorem 3.5 Minimizing the error due to the approximation induced by the nu- 
merical scheme is equivalent to minimizing the norm of the matrices E satisfying 

m- 

Note: Since the linear matricial operator C appears only in the Crank-Nicolson 
scheme, we will restrain our study to the case £ = 0. The generalization to the case 
C ^ can be easily deduced. 

Proposition 3.6 The problem is then the determination of the minimum norm 
solution of: 

MiE + EM2 = F (25) 
which is a specific form of the Sylvester equation: 

AX + XB = C (26) 

where A and B are respectively m by m and n by n matrices, C and X, m hj n 
matrices. 
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3.3 Minimization of the error 
3.3.1 Theory 

Calculation yields: 

The singular values of Mi are the singular values of the block matrix ( 



(27) 

13^ + 5^ (3 [5 + e) 



[3{5 + e) + (3^ 



1. e. 

1 (2/52 + 5^ + s''-{5 + e) ^4/52 + 52 + ^2 _ 25^) (2^ 



of order ^^^^y^, and 



1 (2/?2 + 52 ^ ^2 ^ ^ ^4^2 + 52+^2_25e) (29) 



of order 



The singular values of M2 are a^, of order y? and 7^, of order y- 
Consider the singular value decomposition of the matrices Mi and M2: 

U^M.V, = ( f 1; ) , t/jMir, = ( ) (30) 

where f/i, Vi, ?72, V2, are orthogonal matrices. Mi, M2 are diagonal matrices, the 
diagonal terms of which are respectively the nonzero eigenvalues of the symmetric 
matrices Mi^Mi, Ma^Mg. 

Multiplying respectively on the left side by ^f/i, on the right side by V2, yields: 

MiEV2 + U^ E M2V2 = U^ FV2 (31) 
which can also be taken as: 

^Ui Ml Vi ^Vi E V2 +^ Ui E ^U2 ^U2 M2 V2 = F V2 (32) 



Set: 



V,EV2=\ |ii ^]JU,E^U2=[ ^ ^ I (33) 
F21 F22 ) \ E21 E22 



Fu F12 



We have thus: 



U,FV2=\ tii tii I (34) 
i'21 i'22 



Ml Ell Ml ^12 V / ^11 M2 \ _ I Fii F12 , 

^ 1 I;M2 ) " I ^1 ^2 ' ^''^ 
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It yields: 



MiEu±EiiM2 = F\i 

M^Eu = Fi2 (36) 

E21 M2 = F21 



One easily deduces: 



£"12 = ^ F12 
E21 = F21 M2 



(37) 



The problem is then the determination of the En and En satisfying: 



Ml Eu + Eu M2 = Fu (38) 

Denote respectively by eTj, efj the components of the matrices E, E. 
The problem [38] uncouples into the independent problems: 
minimize ^ 

Cjj + Cij (39) 

under the constraint 

Mui + ^ = Fuij (40) 
This latter problem has the solution: 



Gij — ■ 2 2 



(41) 



The minimum norm solution of [25] will then be obtained when the norm of the 
matrix Fu is minimum. 

In the following, the euclidean norm will be considered. 
Due to 

ll^il 



< lli^ll < \\Ui\\ \\F\\ \\V2\\ < ||f/i|| IIV2II \\MrU,^,act + U,,actM2-Mo\\ (42) 



Ui and V2 being orthogonal matrices, respectively Ux — 1 hy rix — 1, rit hj rit, we 
have: 

||f/if = n,-l , \\V2f = nt (43) 

Also: 

\\M,r = '^^{2(3' + 5^ + e') , ||M2f = |(a2 + f) (44) 

The norm of Mq is obtained thanks to relation f[TU|) . 
This results in: 



li^llll < ^Jnt {Ux - 1) |||f/e.acj (\/^\^ V2/32 + 52+^2+ y^y^^T^) + llMol 

(45) 
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II Full can be minimized through the minimization of the second factor of the right- 
side member of (H5ll . which is function of the scheme parameters. 



||f^exarf|l is a constant. The quantities y ■\/2]3^~+~P~+~e^ , ^foP^^T^ and ||Mo 
being strictly positive, minimizing the second factor of the right-side member of 
can be obtained through the minimization of the following functions: 



UoL^l) = ^MT^ (46) 
f3{a,(3,-f,6,e) = \\Mo\\ 

4 Numerical example: a new DRP scheme 

Consider the scheme (EI) where the values of Px, ^x, and Ex are given by ( |T3l) . 

Let, in a first time, the values of the coefficients a, f3t, 7, 6t, and St remain unknown, 
and advect a sinusoidal signal 

■u = cos I A; (x — ct) ] (47) 

through this scheme, with Dirichlet boundary conditions, (c is taken equal to 1, and 
k = tt). 

Calculation yields then: 



/2(a,7) = _|_ ^2 





1 


2 s 


2 / 




2 










2h 





2 . 1 2 X 2 



(48) 

Minimum values for /i and /s can thus be obtained choosing negative values for /3t, 
while choosing positive ones for St and Et, the absolute values of which are respec- 
tively close to those of Px, Sx and Ex- /2 is minimized choosing 7 = 0. 

In the following, we have choosen to set: 

Pt = -0.9 p^P' 

5t = -0.9 (49) 

^ Et = -0.9 E"J" 

and a = 10. 

The value of the L2 norm of the error, for: 

i. case 1: our new scheme, with cfl = 0.9; 

ii. case 2: the Lax scheme, with cfl = 0.9; 

is displayed in Figure [H The error curve corresponding to the first case is the 
minimal one. 
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Figure 1: Value of the norm of the error. 



5 Conclusion 

The above results open new ways for the building of DRP schemes. It seems that 
the research on this problem has not been performed before as far as our knowledge 
goes. In the near future, we are going to extend the techniques described herein to 
nonlinear schemes, in conjunction with other innovative methods as the Lie group 
theory. 
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